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A model-independent metliod to study tlie possible evolution of dark energy is presented. Optimal 
estimates of the dark energy equation of state w are obtained from current supernovae data from 
[A. G. Riess et al, Astrophys. J. 607, 665 (2004).] following a principal components approach. 
We assess the impact of varying the number of piecewise constant w estimates A^p using a model 
selection method, the Bayesian information criterion, and compare the most favored models with 
some parametrizations commonly used in the literature. Although data seem to prefer a cosmological 
constant, some models are only moderately disfavored by our selection criterion: a constant w, 
w oc a, w oc z and the two-parameter models introduced here. Among these, the models we find by 
optimization are slightly preferred. However, current data do not allow us to draw a conclusion on 
the possible evolution of dark energy. Interestingly, the best flts for all varying-n; models exhibit a 
■u) < — 1 at low redshifts. 

PACS numbers: 98.80.Cq, 95.36.+X 
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I. INTRODUCTION 



The discovery of the late-time accelerated expansion 
of our universe using supernovae (SNe) observations Q 
and independent corroboration coming from additional 
observational tests 2], have often been interpreted as 
indications of the existence of a previously undetected 
dark energy. Furthermore, this appears to be the dom- 
inant constituent of the universe, accounting for 70% of 
its total energy density. Our understanding of the nature 
of dark energy is severely limited by the little we know 
about the dependence on the redshift z of its equation 
of state parameter w{z), defined as the ratio between its 
pressure and its energy density. In practice, observations 
are analyzed while assuming a particular form of w{z)^ 
w = const being the most common one. More elaborate 
parametrizations have been designed to capture the dif- 
ferent behaviors of dark energy models, yet none provides 
a standard description for w. Among their usual deficien- 
cies are: the lack of fiexibility to cope with rapid evolu- 
tion (as pointed out in 0), biasing caused by imposing 
priors over model parameters, and model-dependent re- 
sults. For example 4] warns about the pitfalls of assum- 
ing a constant w or w{z) —1 which, respectively, may 
cause an erroneous reconstruction remarkably consistent 
with observational constraints or an underestimation of 
uncertainties. Studies of possible theoretical motivations 
for w{z)< —1 are available in the literature (see for ex- 
ample p and references therein). 

Principal components (PC) analysis was first used to 
address the problem of parametrizing w{z) in |^], more 
recent examples of the study of dark energy through PC 
being [3, |^ |9|, |10|. In such approach w{z) is described 
in terms of a basis of orthogonal functions whose form 
is dictated by the constraining capabilities of data. A 
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shortcoming of employing PC in this context is that even 
the best determined estimates present broad, difficult-to- 
interpret component functions 0, which obscures the 
intuitive interpretation of the resulting estimates. On 
the other hand, the localized principal components (LPC) 
approach taken by Huterer and Cooray in 7] benefits 
from a more straightforward interpretation: the equa- 
tion of state is now described by uncorrelated and lo- 
calized piecewise constant estimates at different redshift 
regions. Still, the parametrization through LPC is not 
completely model-independent since the redshift bin ar- 
rangement needs to be determined by hand. In Ref. 
the choice was made on the assumption that better con- 
straining capabilities are attained at lower redshifts. As 
will be seen later, we replace this assumption by a de- 
tailed inspection of the possible redshift bin configura- 
tions. Namely, we are asking the following question: Is 
there an optimal structure of the piecewise constant w{z) 
that could facilitate its interpretation, with respect to a 
given dataset? An optimized reconstruction of this sort 
would provide useful information on the phenomenology 
of dark energy, together with a model independent ac- 
count of the constraining power of the dataset under con- 
sideration. 

In this paper we present our efforts to construct such a 
parametrization: In Section^ we show the methodology 
leading to an optimal LPC description of the equation of 
state. First we detail the procedure to obtain a piece- 
wise constant w whose construction is guided by data. 
Since we are unable to decide, on theoretical grounds, 
what number of segments (parameters) to employ in our 
piecewise constant w{z)^ we recur to a model selection 
framework to probe the space of possible parametriza- 
tions. This and other statistical tools used here are also 
described in Section^ Our results and discussion are 
presented in Section Hill where the point is made that 
both the one- and two-parameter optimized LPC mod- 
els are equally suitable descriptions under the model se- 
lection criterion and dataset considered. Concluding re- 
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marks are presented in section HVl 

II. METHODOLOGY 
A. Optimized piecewise constant w{z) 

We choose to represent w{z) by a piecewise constant 
function, assigning a value Wi at each redshift interval 
Zi). This is expressed simply as: 

w{z) = ^Wibi{z), (1) 

where a function bj{z) equals unity within the bin 
(zj-i^Zj) and zero elsewhere. The number of segments 
A^p is defined to be the number of model parameters. 

A model will also be determined by two nuisance pa- 
rameters, the matter density Qrn and the present day nor- 
malized Hubble expansion rate h. In addition, to fully 
determine a model we need to specify the redshift inter- 
val boundaries {^div}- The freedom in choosing the lat- 
ter generates a whole family of parametrizations for fixed 
A^p. For example, when restricted to A^p = 2 we could set 
the boundaries at {0,0.9, 1.8} which would imply differ- 
ent data-fitting capabilities than the choice {0,0.3, 1.8}. 
From now on, z^w refer exclusively to intermediate divi- 
sions since our analyses are all for the (0, 1.8) range. 

In our examples we perform a thorough inspection of 
the fitting properties of all possible {^^div} configurations. 
The optimal model for each Np is identified according to 
the following outline: We start by finding the most likely 
Wi values, with respect to some dataset, for all mod- 
els in a family. We then identify the optimal case, the 
one with highest likelihood (defined below), and refer to 
it as the optimized model for the A^p family, or simply, 
the A/'p-parameter model. This optimization procedure 
fully exploits the descriptive power of piecewise constant 
parametrizations, since it ensures a data-driven recon- 
struction of w. 

The likelihood, C = follows from the statis- 

tic, which is in turn given by 

N / obs _ model/ >)>)2 

where is some observed physical quantity, a its corre- 
sponding uncertainty, y^'^^^^ the quantity calculated as- 
suming some model given by parameters p and the sum 
is over N datapoints. 

B. Localized principal components 

As mentioned in Section ^ our aim is to describe w in 
terms of a basis in which parameter errors are uncorre- 
lated and eigenfunctions are visually easy to interpret. 



The standard PC analysis provides, as an alternative 
to Eq. n an expansion 

^(^) = X]<?iei(z), (3) 

where the new parameters q present uncorrelated errors. 
Nevertheless, the window functions e — the eigenfunc- 
tions of the decorrelation matrix relating the set of orig- 
inal parameters to the uncorrelated ones — are in gen- 
eral difficult to interpret. They are often oscillatory and 
non-zero for the whole redshift range under consideration 
{i.e. not localized). This motivates the use of the LPC 
approach, which yields highly localized and mostly posi- 
tive window functions that allow a direct interpretation 
of the resulting estimates. As pointed out in jl^, this is 
a consequence of choosing the square root of the Fisher 
matrix as the decorrelation matrix. 

We obtain LPC estimates for our piecewise constant 
models following the prescription of J7], as described in 
the Appendix. 

C. Priors 

As mentioned earlier, setting priors on w{z) may de- 
crease our chances of retrieving the correct model. There- 
fore, we only use the top-hat prior implicitly defined by 
the size of the parameter space which, when chosen to be 
large compared to likelihood widths, resembles a 'no pri- 
ors' situation. Later we will see that ranges with consid- 
erable fractions in the w < —1 sector are commonplace. 

Since our ability to constrain w{z) depends strongly 
on our knowledge of ^ , a tight prior will help us 
attain optimal model-recognition capabilities li]. Con- 
sidering this we marginalize using a somewhat optimistic 
fiat prior 0.20 < Qrn < 0.28, in agreement with the la 
constraints of the three-year analysis of WMAP 15]. Ad- 
ditionally, the Hubble parameter is assigned a fiat prior 
0.67 < h < 0.79, compatible with 2a constraints from 
WMAP. A fiat universe is assumed throughout. 

D. Observational test 

The observational quantity used to compute 
through Equation ^ is the one usually reported by ex- 
periments: the distance modulus /i(z) = m{z) — Mobs 
associated to type la supernovae. The apparent magni- 
tude, assumed here to have Gaussian errors, is given by 

m{z) = M - 5 logio ^ + 5 log^o Di^ + 42.38. (4) 

while the Hubble-parameter-free luminosity distance is 
defined as 

^ , ^ du 
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M is the absolute magnitude of supernovae, H{z) 
the Hubble expansion rate and Hq its present day 
value. The perfect degeneracy between M and h = 
10~^Ho km~^ s Mpc is what allows us to reduce the set of 
nuisance parameters to {l^m, h}. An enlightening discus- 
sion on the subtleties involved in pairing the theoretical 
and observed m{z) can be found in Ref. 16]. 

For this analysis we use 156 supernovae from the Gold 
dataset of Riess et a/. 17], which includes 9 SNe at z > 1. 
Our Zdiv arrangements are then restricted to contain at 
least 10 SNe in their last bin, implying Zdw ^ 1- 



E. Model comparison 

As noted in a model selection approach is essen- 
tial when comparing different theoretical models in the 
light of observations; in contrast with just exhibiting es- 
timated values and corresponding confidence levels of the 
model parameters associated to the various models. We 
take this approach to assess the data-fitting merits of dif- 
ferent Np models independently of the number of model 
parameters. 

The problem of choosing the most probable model, 
given certain dataset, finds a simple solution with the 
Bayesian information criterion (BIG, ^H): 

BIG = -21n/:^ax + A^pln7Vdata, (6) 

where >Cmax is the maximized likelihood. Recent exam- 
ples of the use of the BIG in a cosmological context can be 
found in %^,'2d\ amone; other sources. This statistical 
tool was designed as an approximation to the Bayes fac- 
tor B — the ratio of posterior likelihoods used in Bayesian 
statistics to compare a pair of competing models — re- 
sulting in the approximate relation BIG ~ —2\nB. An 
exact solution requires calculating the Bayesian evidence 
|2l|, a more sophisticated task involving integration of 
the likelihood-prior product over the parameter space for 
each model (for an implementation of Bayesian evidence 
to determine the order of a polynomial w{z) see jl4]). 

When priors over model parameters are unknown, as 
is our case, the BIG is a convenient alternative to the 
Bayesian evidence since it is, by design, independent of 
parameter space size. Models may be ranked according to 
the difference ABIG with respect to a base model, chosen 
here to be the 'zero-parameters' AGDM. As listed in Ref. 
I23, a ABIG value between 2-5 (5-10) translates 
into positive (strong) evidence in favor of the model with 
smaller BIG value, while ABIG > 10 implies that the 
evidence is decisive. A ABIG < 2 means neither model 
is preferred by the criterion. 



III. RESULTS AND DISCUSSION 

First let us define AGDM as the base model for model 
selection. The next simplest possibility, constant is 



equivalent to the A^p = 1 model, for which we find an 
estimate wi = —0.90 ± 0.12 (68% confidence level). 

We now illustrate the exploration that identifies the 
optimal bin arrangement for fixed A^p. Fig. Q shows the 
dependence of the bin boundary position z^w for 

A^p = 2 models. We identify the model with z^iv = 0.08 
as the one with a smaller and label it as the 'optimized 
model'. Large-x^ (worse fit) models are limiting cases 
resembling the constant w model, their and preferred 
wi values tend to the ones of the single-parameter case 
(later shown in Table P). This is not surprising if we bear 
in mind that 95% of SNe belong to the first redshift bin 
when Zdiv ~ 1- 
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FIG. 1: Dependence of the bin boundary Zdw that 

characterizes LPC models with two segments. 

In Table ^ we present the models obtained by means 
of optimization for TVp = 1 — 5, together with the base 
model AGDM and two other parametrizations commonly 
used in the literature: w{z) = wq -\- WzZ 23] and w{a) = 
'^0 + '^0(1 — <^) ^24]. We also include the model selection 
quantity ABIG, all statistics are obtained using the Gold 
dataset. 

Model comparison in light of the BIG criterion shows 
that AGDM is the model preferred by data. A collec- 
tion of models that exhibit low but are moderately 
disfavored are Np = 1^ Np = 2^ w oc z^ w oc a and the 
equal errors model explained below. The Np = 3 model 
turns out to be strongly disfavored while A^p = 4 and 
TVp = 5 are decidedly disfavored, all this in terms of the 
rules proposed in [2^ . 

It is worth mentioning that all time- varying w models, 
favored by the BIG or not, suggest a principal component 
value < — 1 at low redshift s. The table also suggests 
that a larger number of LPG estimates implies more ex- 
treme values for Qi. This reminds of certain parametriza- 
tions that, when extended to second order expansions, 
acquire huge values in their expansion coefficients, as 
studied in 0. We won't be concerned by the huge values 
found, since the models in question are already rejected 
by the BIG criterion. Still, it is somewhat disconcerting 
that the further redshift divisions are added, estimates 
are more and more inconsistent with qi ~ —1. 

A more detailed analysis of the optimized Np = 2 case 
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TABLE I: Parameter estimates together with and BIC values for various models. The difference ABIC is calculated with 
respect to the base model ACDM. Redshift arrangements are shown for our optimized models. 



Model 


Optimized {^^div} 
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qi 


q2 




q4. 


q5 Xmin 


BIC 


ABIC 


ACDM 




=-1 












178.6 


178.6 


0.0 


A^p = l 








-0.90 








177.9 


182.9 


4.4 


A^p = 2 


{0.08} 






-2.05 


-0.60 






172.7 


182.8 


4.3 


A^p = 3 


{0.43,0.60} 






-1.08 


2.01 


-1.90 




171.2 


186.4 


7.8 


A^p = 4 


{0.05,0.55,0.70} 






-3.56 


-0.63 


7.45 


-25.01 


170.5 


190.7 


12.1 


iVp = 5 


{0.05,0.50,0.60,0.70} 






-3.62 


-0.67 


9.80 


-3.94 


-40.16 169.6 


194.8 


16.3 


eq. err. 


{0.18} 






-1.28 


-0.46 






173.4 


183.5 


4.9 


It; oc a 




-1.32 


2.61 










173.8 


183.9 


5.3 


W (X z 




-1.27 


1.63 










174.3 


184.4 


5.8 



shows that qi is loosely constrained with respect to q2. 
This is illustrated in Figure [21 where likelihoods (top 
plot) and constraints on model parameters q^ (middle 
plot) are shown. The bottom plot of the figure shows 
the highly localized, mostly positive window functions 
resulting from the LPC approach. 

An alternative A^p = 2 model is obtained by a differ- 
ent optimization criterion which, instead of accepting the 
model with minimum x^, requires errors in parameter es- 
timates to be of comparable size. This case exemplifies 
how different optimization criteria may result in different 
iL'-estimates. In figure |31 we sketch the interplay between 
qi and q2 estimates when 2:div is varied. A detailed anal- 
ysis points to the z^w = 0.18 model as the one with 
matching errors, which is labeled the equal errors model. 
As seen in Table ^ this model has practically the same 
BIC value as the Np = 1 and Np = 2 models. Its cor- 
responding likelihood, parameter estimates and window 
functions are plotted in Figure SI 

It is interesting that two of the models with good fits to 
the data, the optimized Np = 2 and equal errors ones, are 
reminiscent of some more complex models — for example 
the ones in 13, 0, — which suggest w < —1 Sit low 
redshift and w ~ —0.5 elsewere. However, our models 
rely on two parameters only, which makes them more 
attractive in terms of model selection, making a better 
case for the rapidly varying w possibility. 



IV. CONCLUSIONS 

This paper introduces a new framework to study the 
model-independent implications of data regarding the 
possibility of dark energy evolution. Extending the ap- 
proach proposed in 0, we find the optimal set of eigen- 
functions that describe w. We focus on a currently 
available SNe dataset follow a model selection 

scheme, in terms of the Bayesian information criterion, 
which points to ACDM as the preferred model. A class of 
models is found to provide better fits to data but is mod- 
erately disfavored by the BIC. These include our two- 
parameter models: the one that minimizes the 



one with equal errors on both estimates, as well as the 
constant w parametrization. Also in this category, with 
slightly worse BIC values, are the 'linear in redshift' and 
'linear in scale factor' parametrizations. Optimized mod- 
els with more than two parameters are strongly disfa- 
vored. 

In the light of the data considered, mutually exclusive 
phenomenologies are found to be equally favored by the 
BIC, for example, the constant w and optimized A^p = 2 
models which are, respectively, consistent and marginally 
consistent with ACDM at a 95% confidence level. This 
reproduces a problem commonly found in the literature 
when considering flexible forms of w{z): the constant 
w model leads to ~ — 1 at all redshifts, but if w is 
allowed to evolve we find that < — 1 at low redshift 
and w > —1 elsewhere. Our approach may prove useful 
in reconstructing a true model with a rapidly varying 
however, current data give no information on whether 
this option is preferred over a slowly varying or constant 
even when the rapidly varying w is described with a 
simple model such as our two-parameter ones. 

Interpretation of the favored models should be made 
while considering that resulting estimates are "integrated 
properties of w(z)" [2^, and so the underlaying physical 
model may differ substantially from the reconstructed 
one. In particular, the z^iv found in our optimized mod- 
els is not necessarily a proof of change in dark energy 
dynamics at exactly that time, it could be seen as a fea- 
ture of the reconstruction employed. 

Further research, adding other observational tests to 
the analysis, may improve constraints on the equation 
of state. Of special interest are weak lensing measure- 
ments |27|J , or the observed Integrated Sachs- Wolfe effect 
[2^ , which for certain cases provides a probe of the total 
change in w [2^. This is postponed for a future investi- 
gation. 
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FIG. 2: (color online) The optimized A^p = 2 model. Top: 
Likelihoods for the optimized 2-parameter model. Note that 
the low-redshift estimate (red/gray line) is broader than the 
medium- and large-redshift estimate (black line). Color con- 
ventions are the same for the three plots. Middle: Uncorre- 
lated principal component estimates of w{z) (thin, horizontal 
lines) with corresponding 1-cr (soUd) and 2-cr (dotted) rect- 
angles showing parameter constraints (horizontal edges) and 
redshift span (vertical edges). Bottom: Window functions 
relating original piecewise coefficients to LPC ones. 
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APPENDIX 

In this appendix we show, following 7], the procedure 
to decompose w in terms of localized principal compo- 
nents. 

We first calculate the Fisher matrix F (marginal- 
ized over nuisance parameters) for some w model and 
with respect to a given dataset. We then obtain its square 
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FIG. 3: Parameter estimates qi (thick line) and q2 (thin line) 
for varying Zdw and their corresponding 2a confidence regions. 
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FIG. 4: (color online) The equal errors model, a 2-parameter 
model similar to the optimized one but with comparable er- 
rors on both ^-estimates. Likelihoods (top), ^-estimates de- 
scribing w{z) (middle) and window functions (bottom) are 
displayed following the color coding of Figure |21 



6 



root 



(A.1) 



where the diagonahzed matrix of eigenvalues A and the 
orthogonal matrix of eigenvectors O satisfy 



F = O'AO. 



(A.2) 



It is straightforward to check that F^/^ is indeed the 
'square root matrix' of F if we consider the product 

pl/2pl/2^ 

LPC estimates q are given in terms of the original pa- 
rameters p by q = Wp, where W is a normalized version 
of F^/^ for which each row adds up to 1 (Roman indices 
run from 1 to Np and summations are explicitely shown) : 



provide the desired orthogonal basis. Furthermore, when 
considering the covariance of the original parameters 
(ApiApj) = and the diagonality of W, the new 

errors Aqi are found to be uncorrelated, their covariance 
being 



{AqiAqj) = 



k I 



(A.4) 



Our expansion in terms of the window functions is then 
written as: 



(A.5) 



F 



1/2 



(A.3) 



The rows of W then define the window functions e which 



where instead of ordering the principal components, as 
is usual, according to the magnitude of individual uncer- 
tainties, we respect the original indexing to maintain the 
transparency of the interpretation. 
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